## moving window.

####Load the libraries

 #list the libraries needed
necessary=c("seewave",'sound', "tuneR", "geoR", "caTools", "sp", "rgdal", "adehabitat")
?adehabitat
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)

#load the libraries
for (lib in necessary) library(lib,character.only=T)
    
##set the working directory

setwd ("C:/BIOACOUSTICS/movingwindow/")

outdir = ("C:/BIOACOUSTICS/movingwindow/")

### load in the sounds  ###########

library(tuneR)
s6 <- readWave("oriole_small.wav")
play (s6)
s6

s7= loadSample("oriole_freq_range.wav", filecheck=TRUE)

specASCII = spectro(s7,f=16000, plot = F)   ## do Fast Fourier Transfer on .wav sample and punch out a matrix
specASCII [1:3,]   # have a look at it

t1 = as.matrix (specASCII)  # turn the spectogram into a matrix
t2 = t(t1)   #transpose the matrix to view the data as we see it in the spectogram.
t3 = as.asc (t2)   # turn the transposed matrix into an ascii file
x = (1:nrow (t3))
y = (1:ncol(t3))
image (t3, xlim = range(x, finite = TRUE), ylim = c(30, 80), col = terrain.colors (20))    #look at an image of the ascii to see where the values sit
locator (1)
locator (1)

## some notes on analysis:  we need to convert the numbers into a dimensionless form ie like transforming numbers for a similarity matrix.  Each spectogram needs
# to represent not the unique numbers (that will be based on how close the animal is to the microphone) but their relative dimension compared to other species.  
#  i.e. the 3 dimensional shape of the sound will remain relatively similar just change it's size based on the strength of the call.  If all calls are reduced
# to their relative shape then the shapes can be compared rather than the decibels.


# playing with moving windows without knowing where to look for the real stuff.
xl = c(15, 50)
yl = c(50, 66)

t4 = subsetmap(t3, xlim = xl, ylim = yl)  ## first window
ty =as.matrix (t4)
write.csv (ty, "ty.csv")

image (t4, xlim = xl, ylim = yl, col = terrain.colors (20)) # have a look

t5 = subsetmap (t3, xlim = xl + 10, ylim = yl) ## the next chunck of map

hist (t5, density = 20, angle = 70)

image (t5, xlim = xl + 10, ylim = yl, col = terrain.colors (20))    # have another look 

t6 = subsetmap (t3, xlim = xl + 30, ylim = yl) ## and now the top 
image (t6, xlim = xl + 30, ylim = yl, col = terrain.colors (20)) # how does it compare.

t7 = subsetmap (t3, xlim = xl + 40, ylim = yl) ## and now the top 
image (t7, xlim = xl + 40, ylim = yl, col = terrain.colors (20)) # how does it compare.

t8 = subsetmap (t3, xlim = xl + 50, ylim = yl) ## and now the top 
image (t8, xlim = xl + 50, ylim = yl, col = terrain.colors (20)) # how does it compare.

t9 = subsetmap (t3, xlim = xl + 60, ylim = yl) ## and now the top 
image (t9, xlim = xl + 60, ylim = yl, col = terrain.colors (20)) # how does it compare.

t10 = subsetmap (t3, xlim = xl + 190, ylim = yl) ## and now the top 

image (t10, xlim = xl + 190, ylim = yl, col = terrain.colors (20)) # how does it compare.



t11 = subsetmap (t3, xlim = xl + 200, ylim = yl) ## and now the top 
image (t11, xlim = xl + 200, ylim = yl, col = terrain.colors (20)) # how does it compare.

t12 = subsetmap (t3, xlim = xl + 210, ylim = yl) ## and now the top 
image (t12, xlim = xl + 210, ylim = yl, col = terrain.colors (20)) # how does it compare.

t13 = subsetmap (t3, xlim = xl + 505, ylim = yl) ## and now the top 
image (t13, xlim = xl + 505, ylim = yl, col = terrain.colors (20)) # how does it compare.

t14 = subsetmap (t3, xlim = xl + 510, ylim = yl) ## and now the top 
image (t14, xlim = xl + 510, ylim = yl, col = terrain.colors (20)) # how does it compare.

t15 = subsetmap (t3, xlim = xl + 520, ylim = yl) ## and now the top 
image (t15, xlim = xl + 520, ylim = yl, col = terrain.colors (20)) # how does it compare.


hist (t4, col = "black", border = "white")
hist (t10, density = 15, angle = 120, col = "yellow", add = T)
hist (t14, density = 20, angle = 70, col = "grey", add = T)
hist (t7, density = 5, angle = 20, col = "white", add = T)



tt = "mean" = NULL
tt$meant4 = mean (t4)
tt$meant5 = mean (t5)
tt$meant6 = mean (t6)
tt$meant7 = mean (t7)
tt$meant8 = mean (t8)
tt$meant9 = mean (t9)
tt$meant10 = mean (t10)
tt$meant11 = mean (t11)
tt$meant12 = mean (t12)
tt$meant13 = mean (t13)
tt$meant14 = mean (t14)
tt$meant15 = mean (t15)

ttt = as.data.frame (tt)

?abline
tsd = NULL

tsd$sdt4 = sd (t4)
tsd$sdt5 = sd (t5)
tsd$sdt6 = sd (t6)
tsd$sdt7 = sd (t7)
tsd$sdt8 = sd (t8)
tsd$sdt9 = sd (t9)
tsd$sdt10 = sd (t10)
tsd$sdt11 = sd (t11)
tsd$sdt12 = sd (t12)
tsd$sdt13 = sd (t13)
tsd$sdt14 = sd (t14)
tsd$sdt15 = sd (t15)

tsdData = as.data.frame (tsd)

edit (tsdData)




#get some info on each window
mean(t4)
mean(t5)
mean(t6)


########### create a filled contour of the area of interest to see how things work.

x = (1:nrow (t3))
y = (1:ncol(t3))
z = t3

filled.contour(x = (1:nrow (t3)),
y = (1:ncol(t3)),
         z = t3,
         nlevels = 20,
         labels = NULL,
         xlim = range(x, finite = TRUE),
         ylim = c(30, 80),
         zlim = range(z, finite = TRUE),
         drawlabels = TRUE, method = "flattest",
         col = terrain.colors(19))
         
locator (2)     ### locate the coordinates that bound the upper left, upper right, lower right and lower left (create a box for the moving window).


## set up moving window analysis of FFTransformed data

x = t3
k = 14

t4  = runsd(x, k, center = runmean(x,k), 
                 endrule=c("sd"))   


### have a look at the structure of the wav in many different interesting ways

spectro3D(s7,f=16000,wl=500,ovlp=75,zp=50,maga=50,palette=rev.terrain.colors)
ama(s7,f=16000,wl=1024,identify=TRUE)

spectro

# simple plots

spectro(s7,f=16000,osc=TRUE)
spectro(s7,f=16000,scale=FALSE)
spectro(s7,f=16000,osc=TRUE,scale=FALSE)
# manipulating wl
op<-par(mfrow=c(2,2))
spectro(s7,f=16000,wl=256,scale=FALSE)
title("wl = 256")
spectro(s7,f=16000,wl=512,scale=FALSE)
title("wl = 512")
spectro(s7,f=16000,wl=1024,scale=FALSE)
title("wl = 1024")
spectro(s7,f=16000,wl=4096,scale=FALSE)
title("wl = 4096")
par(op)
# vertical zoom using flim
spectro(s7,f=16000, ylim=c(2,6))
spectro(s7,f=16000, ylimd=c(2,6))
csh(s7,f=16000,wl=512,ovlp=50)
dynspec(s7,f=16000,wl=1024,ovlp=50,osc=TRUE)
dev.off()

timer(s7,f=16000,threshold=5,smooth=50)
   
# get the spec entropy for the sound sample and output it as a csv for comparison  
           
test=csh(s7, f=16000, wl = 82, wn = "hanning", ovlp = 50, threshold = NULL,
plot = FALSE)
 
 as.data.frame(s7, row.names = NULL, optional = FALSE,  
              stringsAsFactors = default.stringsAsFactors())


